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ABSTRACT 



A numerical study was performed to investigate the dynamic 
crack propagation in fibrous composite plates utilizing the 
finite element method. A rectangular plate of uniform 
thickness, which had a propagating central crack, was used for 
the study. The plate was a unidirectional composite panel and 
the load was applied in the longitudinal direction of the 
composite plate. Fracture energies were calculated for given 
speeds of cracks. The objective of the study was to examine 
the effect of different composite material properties, crack 
speeds, and densities on the fracture energy. The 
discontinuous node- release technique was used to model the 
crack propagation. 

The numerical study showed that the fracture energy was 
higher for a lower elastic modulus ratio, if all other 
conditions were held the same. Furthermore, a lower crack 
propagation velocity or a lower material density resulted in 
a higher fracture energy, respectively, provided the rest of 
the parameters held constant. 
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I. 



INTRODUCTION 



Dynamic fracture deals with problems having a crack 
within a body where the inertial forces play an important 
role. The early investigation of the dynamic fracture analysis 
goes back to the case of unstable crack propagations in some 
ships during World War II. The dynamic fracture problems 
involve crack initiation and crack propagation. There are two 
principal types of problems. One is applying a dynamic load to 
a body containing stationary cracks, and the other is 
concerned about a body with a moving crack. The second type of 
problems can be classified into two different modes; 
propagation and generation modes. 

Stationary crack problems involve cracks within structures 
that are subjected to dynamic loading, such as an impact load. 
For example, Chen performed a numerical analysis to find the 
dynamic stress intensity factor for a rectangular plate having 
a central crack. He used the finite difference technique and 
discussed the dynamic behavior of the stress intensity factor. 
He argued that the oscillations on the stress intensity versus 
time curve were mainly due to scattering phenomena from the 
crack tip and the boundary surfaces. He concluded that the 
overshooting of the dynamic stress intensity factor is 
attributed mainly to the geometry and loading. [Ref. 1] 
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The generation mode involves moving cracks expanding at 
prescribed rates. The main objective of the study is to 
compute the stress intensity factor or the fracture energy and 
to determine the fracture toughness. The propagation mode 
deals with the cases where fracture toughness is known. The 
objective is to predict the crack propagation and its arrest 
for a given load. One of useful models to investigate the 
crack propagation and arrest was the double cantilever beam. 
Kanninen and his coworkers [Ref. 2] proposed the model and 
studied it extensively. Some of studies in the dynamic 
fracture were listed in the references [Ref. 3-7] . 

Although there are many literatures available for the 
dynamic fracture of isotropic materials, the studies of 
composite materials are not so common. Some of the fracture 
studies of composites were shown in references [Ref. 8 -10]. 

This study considered a rectangular plate made of 
unidirectional composite and containing a central, moving 
crack subjected to a uniform tensile load. The finite element 
method utilizing bilinear rectangular elements was used for 
this study. A parametric study to compare fracture energies in 
composite plates for various conditions was conducted. 
Different elastic modulus ratios, crack velocities, and 
material densities were considered. 
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II. 



MATHEMATICAL DERIVATION 



A. GOVERNING EQUATIONS 

The governing equations for the problem of elasticity are 
equations of equilibrium, constitutive equations, and 
kinematic equations. For the two-dimensional problem, they are 
given below. Equations of equilibrium are 



where a x , a y , and are the two normal stresses in the x and 
y direction and the shear stress, respectively, p is the 
material density, and u and v are the displacements in the x 
and y direction, respectively. In addition, x and y are 
rectangular coordinate axes and t denotes time. 

Constitutive equations are also known as stress -strain 
relations. The relations are given in a matrix form like 




dx dy K at 2 



( 1 ) 



{<j}= [D] (e) 



( 2 ) 



are, respectively, stress and strain 



where {a} and {e} 
vectors, that is 



{<j}={<j x o y xJ T 



( 3 ) 



and [D] is the material property matrix of size 3x3. For an 
isotropic material, the [D] matrices for plane stress and 
plane strain conditions, respectively, are 



[D] =■ 



(1-u) 2 



1 i) 0 

u 1 0 

1-u 



0 0 



(4) 



and 



[D] =■ 



Ed-u) 



(l+u) (l-2u) 



1 

u 



u 



1-u 

1 



1-u 
0 0 



0 

0 

l-2u 



2 (1-u) 



(5) 



For a composite material, [D] matrix becomes 



[D] = 















^22 


^23 




( 6 ) 


D 21 




^ 33 . 







The details of Eq. (6) are shown in references. [Ref. 11] 

Kinematic equations are also known as strain-displacement 
relations. The infinitesimal strains are expressed in terms 
of displacements as shown below: 
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B . BOUNDARY CONDITIONS 



On the essential boundaries displacements u and/or v are 
prescribed, and on the nonessential or natural boundaries 
tractions are prescribed. Usually on a part of the boundary 
tractions are prescribed, while on the remainder of the 
boundary displacements are known. The tractions are given by 



$*=<W-W 2 y 



(8) 



where <£ x and <£ y are the tractions in x and y directions, 
respectively, and n x and riy are, respectively, the components 
of outward unit normal vector to the boundary. 

Applying the method of weighted residual technique to 
the equations of equilibrium yields 



j 2 =Jq(-P 



d 2 u + 


da x 
A + . 


foxy 


at 2 


dx 


By 


d 2 V 


& X xy + 


d(J y 


dt 2 


dx 


dy 



( 9 ) 



where and w 2 are test functions and Q is the domain of a 
given problem. Applying the weak formulation results in 



-Jc 



d 2 U 

P <*>1 

K 8t 2 1 dx 



0G) 1 

^ v^dO + f - 



0 2 v 

P- — :«,+T 



00 ), 00 ), 

— - +a 

w dx y dy 



0 ( 10 ) 



dt 2 2 

and r is the boundary of domain Q. Rewriting Eq. (10) gives 



Jq(P 



0)! 0 

0 o). 



0 2 U 



dt 2 
02 V 



0t ‘ 



> + 



00) 1 
dx 



0 



00 ). 



00) x 



dy 

00 ), 



0y 0x 



M> dQ - /rfcvdr (11) 



’*yj 



<"4< 

,“> 2 j 



where 



Jq 







02 u 




JqP 


0) x 0 
0 0) 2 


dt 2 

02 V 


>dQ 






.dt 2 




' 0O) t 

dx 


0 


d(0 1 

dy 


<V 

rt 


0 


0O) 2 

dy 


00) 2 
<9x 


y ( 

W 



VdQ 



( 12 ) 



( 13 ) 



and 




( 14 ) 



are the mass matrix, stiffness matrix, and load vector, 
respectively. 
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A finite element discretization is applied to the domain 
and the following matrices are computed at the element domain. 

C . MASS MATRIX 

The mass matrix is given by 



in* P 



(0 



0 



1 

0 Cl) 



dt 2 

dt 2 



}dQ 



( 15 ) 



In this study both bilinear rectangular and linear triangular 
elements will be used to produce the finite element mesh. 
First the finite element derivation will be shown in detail 
using the bilinear rectangular element. The rectangular 
element is shown in Figure 1. 




Figure 1 Bilinear Rectangular Element 
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The shape functions are 



tf 1 = -TT— (b-x) ( c-y ) 
1 Abe 

= (b+x) (c-y) 

z Abe 

tf 3 = — (b+x) (c+y) 
Abe 

# 4 = — (b-x) (c+y) 
Abe 



( 16 ) 



where b and c are greater than zero. The displacements in the 
x and y directions, respectively are interpolated as shown 
below 



u=H 1 u 1 +H 2 u 2 +H 2 u 2 +H i u i 
v=H 2 v x +H 2 v 2 +H 2 v 2 +H i v i 



( 17 ) 



where 



u i 

U 2 



“3 

^3 

U 4 

^4 



is the displacement vector at the nodes of the element . 

Because the shape functions are independent of time the 
acceleration terms become as shown below: 
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c^u 

dt 2 

d^v 



dt 2 



= K fl i 



+// 2 u 2 +// 3 u 3 +// 4 uJ 



0 H 2 0 H 3 0 H i 0 
0 H ± 0 H 2 0 H 3 0 H a 



u. 



u. 



U-. 



U A 



( 19 ) 



Substituting Eq. (19) back into Eq. (15) and applying 
Galerkin's method results in 



Jq s P 



*1 


0 


0 






0 


0 




*3 


0 


0 


*2 


*4 


0 


0 


H, 



H x 0 H 2 0 H 3 0 H a 0 
OH, 0 H 0 0 0 H. 



dCi 



u. 



u. 



u. 



U £ 



( 20 ) 



Finally Eq. (20) becomes 
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Jq«P 



Hi 0 H x H 2 0 H 3 H 3 o h 1 h a o 

0 Hi 0 H x H 2 0 H X H . 3 0 H X H A 

H x H 2 0 Hi o h 2 h 2 o h 2 h a 0 

0 Hj'H . 2 0 H 2 0 H 2 H 2 0 H 2 H 4 

h x h 2 0 H 2 H 3 0 h 3 o h 3 h a 0 

0 Hj_H 3 0 H 2 H 3 0 Hi o h 3 h a 

h x h a o h 2 h a 0 h 3 h a 0 H a 0 

0 h a h a 0 h 2 h a 0 H 3 H a 0 Hi 




(21) 



Carrying out the integrations gives symmetric consistent 
mass matrix as shown below: 



where t is the thickness of the element. Equation (22) is 
called the consistent mass matrix. 

Another choice for the mass matrix is the lumped mass 
matrix. The lumped mass matrix is a diagonal matrix so that it 
reduces computation time and space needed for storage. The 



[M c \ =-^bct 



40201020 
4 0 2 0 1 0 2 

4 0 2 0 1 0 

4 0 2 0 1 

4 0 2 0 

4 0 2 
4 0 
4 



( 22 ) 
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lumped mass matrix is 



10 0 
1 0 
1 



[Mj] =pbct 



0 

0 

0 

1 



0 0 0 0 
0 0 0 0 
0 0 0 0 
0 0 0 0 
10 0 0 
10 0 
1 0 
1 



( 23 ) 



D. STIFFNESS MATRIX 

The stiffness matrix is given by 



io 8 



do) 1 
dx 



doi 1 



do. 



dy 

3o) 0 



dy dx _ 



d Q 



( 24 ) 



Substituting constitutive equations, Eq. (2) , into Eq. (24) 
yields 



do^ 




0 



0 



3<0 2 

dy 



d 

~dy 

dx 



1 


[ e *| 


[£>]j 




1 


\fxy\ 



dQ 



( 25 ) 
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Using kinematic equations, Eq. (25) becomes 



Jo® 



' d(o 1 


0 


aw.^ 




au 

dx 


dx 


ay 


(£>]- 


dv 




a w. 


aw. 


dy ' 


0 


z 

ay 


z 

dx a 




av + au 
ax ay 



dQ, 



( 26 ) 



The kinematic equation can be written in terms of the nodal 
displacements as below: 



du 




dx 




dv 




dy 


► = 


dv + du 




dx dy 


L 



± o -A 

dx dy 



o -A JL 

dy dx 



‘_3_ 

ax 


0 


_a_' 

ay 


o 

aT 

o 


o 

ar 

o 

ar 


n 


_a_ 


_a_ 


af 

o 

a? 

O 


a? 

o 

o 


u 


ay 


dx 







dH x 

dx 



dx 



dH ± 



0 

dH , 






a//. 



dHj 

dx 

0 



a//. 



6y dy 

dH x dH ± dH 2 dH 2 dH 3 dH 3 dH 4 

dy dx dy dx dy dx dy dx 



dy 

dH. 



u. 



u„ 



u. 



u A 



(27) 
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where the last matrix is called [B] matrix and the last vector 
called the displacement vector {u}. Substituting back Eq. (27) 
into Eq. (25) yields 



Jq« 



3o 1 

dx 

0 



3(0. 

0 “5^ 

dy 

3g> 2 3g ) 2 



dy dx _ 



[D] [ B ] dQ{u) 



Using Galerkin's method the first matrix in Eq. (28) 



d(ji> 1 

dx 

0 



0 

3 ( 0 . 



3w 1 



dy 

3(o, 



dy dx 



dH 1 


o 


3// x 


dx 




3y 


0 


dH ± 


dH t 




dy 


dx 


dH 2 


0 


dH 2 


dx 




dy 


0 


dH 2 


dH 2 


dy 


dx 


dH 3 


0 


dH 3 


dx 




dy 


0 


dH 3 


dH 3 


3y 


dx 


dH, 


0 


dH, 


dx 


dy 


0 


3// 4 


dH, 




3y 


dx 



which is equal to the transpose of B matrix, 



( 28 ) 



becomes 



( 29 ) 
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Finally, the stiffness matrix is 



[K] =f QS [B] t [D] [B]dQ. (30) 

Computing the [K] matrix using bilinear rectangular element 
gives 



f-d-blB] t [D] [ B ] tdxdy 



(31) 



where t is the thickness, and 



IB] 



4bc 



( c-y) 0 ( c-y) 0 ( c+y) 0 - ( c+y) 0 

0 -( b-x ) 0 -( b+x ) 0 ( b+x ) 0 ( b-x ) 

.-( b-x ) -(c-y) -(b+x) (c-y) (b+x) (c+y) (b-x) (c+y) 



(32) 



For a general material property matrix [D] of Eq. (6) , the 
element stiffness matrix becomes 



m = 



*11 


*12 


*13 


*14 


in 


*16 


*17 


*18 


*21 


■^22 


*23 


*24 


*25 


*26 


*27 


*28 


*31 


*32 


*33 


*34 


*35 


*36 


*37 


*38 


*41 


*42 


*43 


*44 


*45 


*46 


*47 


*48 


*51 


*52 


*53 


*54 


*55 


*56 


*57 


*58 


*61 


*62 


*63 


*64 


*65 


*66 


*67 


*68 


*71 


*72 


*73 


*74 


*75 


*76 


*77 


*78 


*81 


*82 


*83 


*84 


*85 


*86 


*87 


*88 



(33) 



The detail of K^j is shown in the Appendix. 
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If a material is isotropic and the plane strain condition 



applies, the stiffness matrix is 

[K] = [K n ] + [K s ] . 



[KjJ is given by 



b 3 c 4 ab 2 c 2 y£> 3 c 


-4 ab 2 c 2 


-y J b 3 C 


-4ab 2 c 2 


-J 6 ^c 


4 ab 2 c 2 


-y-jbc 3 4ajb 2 c 2 


-fbc* 


-4ab 2 c 2 


-y-j be 3 


-4ajb 2 c 2 


yjbC 3 


16 

-y £ J C 


-4 ab 2 c 2 


-16 

— 3“^ c 


-4ajb 2 c 2 


-y-jb 3 c 


4a£ 2 c 2 




Ifbc' 


4 ab 2 c 2 


yj be 3 


4ai) 2 c 2 








Ifb'c 


4ai? 2 c 2 


yib 3 C 


-4<3jb 2 C 2 










4a£ 2 c 2 


-l 6 he* 










ifb 3 c 


-4ab 2 c 2 






where 



ff(l-u ) t 

(4bc) 2 (1+u) (l-2u) 



(1-u) 



and [K s ] is given by 



( 34 ) 



( 35 ) 



( 36 ) 
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^-b 3 c 4b z c 2 


-|i> 3 c 


-4b 2 c 2 


-b 3 c 


-4b 2 c 2 


-1 6 £> 3 c 


4b 2 c 2 


±*bc 3 

3 


4jb 2 c 2 


~ 16 bc 3 


-Ab z c 2 


^-bc 3 


-4jb 2 c 2 


-|i>c 3 




16 p. 

-^-b*c 


-4b 2 c 2 


” 16 b 3 c 


-4b z c 2 


^-b 3 c 


4 i> 2 c 2 






lj-bc 3 


4b 2 c 2 


-|i?c 3 


4b z c 2 


^Y~bc 3 








4pb 3 c 


4b 2 c 2 


~b 3 c 


-4b 2 c 2 












4 jb 2 c 2 


~ 16 bc 3 



■^■b 3 c -4£> 2 c 2 

M . be 3 



where 




Et 

(4bc) 2 ( 1 +i) ) 



( 38 ) 



E . LOAD VECTOR 

The load vector is given by 




Linear shape functions can be used to find nodal load values 
due to the traction. This results in lumping tractions into 
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two nodal points associated with the boundary of the element. 
If there is no traction, the load vector is zero at the 
associated nodal points. Linear shape functions are 



where s is the coordinate axis along the boundary, and hi is 
the length of the boundary of the element. Combining all of 
the above procedure gives an equation of matrix, 



which describes the equation of motion with no damping. 
Damping was neglected in Eq. (41) . 

F. LINEAR TRIANGULAR ELEMENT 

Since this study deals also with the linear triangular 
element, the derivation of the matrices will be shown here, 
but not in detail since the detail was shown before in the 
bilinear rectangular element case. The linear triangular 
element, as shown in Figure 2, has three nodal points which 
have coordinate values ( x i»yi) / ( X 2'Y2) > and ( X 3 'Y 3 ) 
respectively. Linear triangular element has two degrees of 
freedom per node. 



JL 

^i =s i*l~ s i 




( 40 ) 



[M] (ii}+ [Jfl {u)=Ir} 



( 41 ) 
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Figure 2 Linear triangular element 



Note that three nodal points are numbered in the counter- 
clockwise direction from an arbitrarily selected corner. 

The shape functions are 



Hi= Ya [a i +b i x+c i y ] 

H 2 = YA [A 2 +B 2 X+C 2Y ] (42) 

tf3=H IA,+B,x+C 2 y] 

where A is the area of the triangle, in other words 



and 




1 Yi 
1 *2 
1 x 3 Y 3. 



(43) 
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( 44 ) 



=*2-^3 "*3^2 B l=y 2 -y3 

^ 2 = * 3 yi -^y 3 b 2 = y 3 -yi 

^3 =^2-^1 s 3 = yi-y 2 



C l=^ 3-^2 
C 2 =*l -*3 
C 3 =x 2 -Xi . 



It can be noted that 



A={A 1 +A 2 +A Z ) /2 . 



( 45 ) 



The displacements in the x and y directions are given by 



u=// 1 u 1 +i/ 2 u 2 +i/ 3 u 3 

v=H 1 v 1 +H 2 v 2 +H 2 v 2 



( 46 ) 



Carrying out similar calculations and procedure gives the 
symmetric consistent, and lumped mass matrices, respectively 



2 




0 10 10 
2 0 10 1 
2 0 10 
2 0 1 
2 0 
2 



( 47 ) 



1 0 
1 




0 

0 

1 



0 0 0 
0 0 0 
0 0 0 
10 0 
1 0 
1 



( 48 ) 



The stiffness matrix is obtained from a combination of two 
components 



19 





q qqu 


b 1 b 2 


fliC 2 u 


qq 


qqu 




cl 


qqu 


Ci^2 


qqu 


qq 


EV 




q 2 




qq 


q qi> 


4 A 2 (1-U 2 ) 






cl 


qqu 


qq 










si 


qqu 












c 3 2 



( 49 ) 



and 



c i 



\K 1 g v 

S 8A 2 (1+u) 



qq qq qq qq qq 
q qq qq qq qq 

^2 q *®2 ^ 2^3 ^ 2-®3 

q 2 qq qq 
q 2 qq 
s 3 2 



where V is the volume of the element, that is, 



V-At 



( 50 ) 



( 51 ) 



and t is the thickness of the element. 

G. TRANSIENT ANALYSIS 

In this study the central difference method is used as a 
time integration scheme. The equation of motion at time t is 
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[M\ {u) c + [C] (u} c + [K] {u) c =tR) c . 



( 52 ) 



The acceleration at time t is given by 

{u} c =— — [{u} c ' aC -2{u} c +{u) c+aC ] . 

At 2 

The velocity at time t is 



( 53 ) 



(u}t=_JL_ [-{u}c-At + { u }t + . C] . (54) 

2a t 

The error associated with the acceleration and velocity 
calculations are of order (At 2 ) . Substituting the relations 
for acceleration and velocity into the equation of motion. 



(At) 



[M] +-i- [C] ){u) t+i£ = 



2 A t 



(l?} t -( [ic] - 



(At) 



[M] )iu) t -(- 



(At) 



[M] --i- [C])W 
2At 



( 55 ) 



t-At 



The method requires the solution at time -At. Using the 
equations of acceleration and velocity, letting time to be 
equal to zero and eliminating {u} At yield 



(u}' aC ={u}° -a t{u)° + — (u}° 

2 



( 56 ) 



Note that {u}° can be found using the equation of motion, Eq. 
(52) , at time zero since {u}° and {u}° are known. 

The central difference scheme is conditionally stable and 
requires that the time step At should be smaller than a 
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critical time step At cr . The critical time step is given by 



where T min is the smallest period. 

Therefore the method can be summarized as; 

1. Compute mass, damping, and stiffness matrices. 

2. Compute {u}° from 



[M] {u}°=Ir}°- [C] (ti}° + [K] (u}° 



where {u}°,and {u}° are initial conditions. 

3 . Compute 

{u}' lt =(ii 0 -At{ii) 0 +- . 

2 

4. Compute effective mass matrix from, 



[M] =■ 



(At) 



[M] 



2At 



[C] 



5. Compute the effective load vector at time t. 



{£} c =te} c -( [K] - — [M] )(u} c -( — ^ - [M] [C] 

(At) 2 (At) 2 2At 



6. Solve for displacements at time t+At from, 

[M] {u} c+aC ={j?} c . 



7. Compute accelerations and velocities at time t, 



(57) 



(58) 



(59) 



(60) 



(61) 



( 62 ) 
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1 



{ u ) c = — — [{u} c aC -2{u) c +{u} c * At ] 
At 2 

{u}n = _l_ [.{u}e-« + {u}f*t] 

2 A C 



( 63 ) 



Note that steps 5-7 have to be performed repeatedly for each 
time increment . 



H. DYNAMIC FRACTURE ANALYSIS 

There are many parameters to describe the crack tip 
behavior for determination of the crack propagation and 
arrest. One of the parameters is the stress intensity factor, 
K, at a crack tip which is defined as 



where 



K= 1 im r _ >oV ^2 nra n (r, 0, t) 



(64) 



°n~ a n ( r ' 6' t) (65) 

is the stress component, and is evaluated near the crack tip. 
In this case, crack propagation is governed by stress 
intensity factor K, and the material property K c , called the 
fracture toughness. A crack becomes unstable if the stress 
intensity factor reaches the fracture toughness of the 
material . 



23 



The dynamic analysis requires K be considered in two 



parts. For crack initiation, 

K(a, 0 't)=K D ( 6 ) (66) 

where a is the crack length, a is the applied stress, t is 
time, and o represents the loading rate. For a propagating 
crack, 



K(a,o, t) =K d (a) (67) 

where a indicates the crack speed. [Ref. 7] 

An alternative criterion for crack motion is, 

G(a,o , t) =R(a) (68) 

where G is the dynamic energy release rate and R is the energy 
dissipation rate required for crack growth which is the 
critical value of G. 

The driving force or the dynamic energy release rate for 
crack propagation has contributions from strain energy, 
kinetic energy, and work done on the body. The driving force 
is the net change in these components per unit area of crack 
extension. The equation for G is, 
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( 69 ) 



r _l , dW_ dU_ dT ) 
b da da da 
1 ( dW_dU_dT) 
bd dt dt dt 

where U is the strain energy, T is the kinetic energy, W is 
the work done by the external loading, a is the crack length, 
and b is the plate thickness at the crack tip. 

The dynamic energy release rate can be directly connected 
to the dynamic stress intensity factor. For plane strain 
conditions, this relation is 

G=±^l A (a)K 2 ( 70 ) 

E 

where E and v, as usual, are the elastic modulus and Poisson's 
ratio, respectively, and A is a geometry independent function 
of the crack speed given by 



(^r) 2 {l~— 2 ) 2 

C 2 C 2 



(1-V) [4(1-^) 2 (1-^) 2 -(2-l^) 2 ] 

c 2 cl cl 



( 71 ) 



where C x and C 2 , respectively, are longitudinal and shear 
(transverse) wave speeds in the material. Wave speeds are 
given by, 




(iC b +4G s /3) 

P 




( 72 ) 
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and 



C 2 ={G s / p ) 1 ' 2 



( 73 ) 



where K b and G s are bulk and shear modulus, respectively. 
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III. 



RESULTS AND DISCUSSION 



A. VERIFICATION PROBLEMS 

In order to verify the developed computer code, two 
example problems with known solutions were analyzed for 
isotropic plates. One was a stationary crack problem subjected 
to a rapidly varying dynamic load, and other was a propagating 
crack problem with a constant speed and subjected to a uniform 
load. 

The first verification problem was a stationary crack 
problem known as Chen's problem [Ref.l]. A rectangular bar 
with a centrally located crack, shown in Figure 3, was 
considered. A uniform tension of magnitude 0.004 Mbar with 
Heaviside- function time dependence was applied. The material 
properties were given as below; shear modulus of 0.76923 Mbar, 
bulk modulus of 1.66667 Mbar, Poisson's ratio of 0.3, and 
material density of 5 gr/cm 3 . Chen solved the problem using 
the finite difference technique. 

Both Chen's solution and the present solution are shown in 
Figure 4. A good agreement between the two solutions was 
obtained. 

As the second example a rectangular plate with a central, 
moving crack was considered [Ref. 3] . The crack was assumed to 
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P(t) 



Figure 3 The Geometry of the Stationary Crack Problem 
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CO 




Figure 4 The Results of the Stationary Crack Problem 
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propagate at a constant speed with uniform stress field of 
68.9 MPa. The crack speed was assumed to be 33% of the 
dilatational wave speed in a steel which is 5904 m/s. The 
material properties were elastic modulus of 206.7 GPa, 
Poisson's ratio of 0.30, and material density of 8000 kg/m 3 . 
The geometry of the problem is shown in Figure 5. [Ref. 3] 

A good agreement was obtained between the analytical and 
the present solutions. Both the analytical and present 
numerical analysis results are shown in Figure 6. 

B. RESULTS OF COMPOSITE PLATES 

A composite plate with a central, moving crack, with a 
similar configuration as the previous example, was 

investigated for different cases. A unidirectional composite 
plate was considered. 

For all cases calculations were made for work done on the 
system, internal energy, kinetic energy and fracture energy. 
The objective of this study was to compare the fracture 
energies for various conditions. The first case was the 
variation of the ratio of elastic moduli of the longitudinal 
and the transverse directions while keeping the others fixed. 
The next case changed crack velocities and, the last case 
varied material densities while keeping the other parameters 
fixed. The geometric and material properties of the problem 
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Figure 5 Finite Element Mesh for the Moving Crack Problem 
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Figure 6 The Results of the Moving Crack Problem 
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are length of 3m, width of lm, elastic modulus ratio in the 
transverse direction of 10 GPa, Poisson's ratio in the 
transverse direction associated with loading in longitudinal 
direction of 0.25, and uniform loading in the longitudinal 
direction of 70 MPa. 

Figure 7 shows work and energy variations during a crack 
propagation. For all cases throughout this study, work was 
higher than internal and kinetic energies, and internal energy 
was higher than kinetic energy. 

In the first case, the elastic modulus ratios were 
changed. Three different ratios between the elastic moduli in 
the longitudinal and transverse directions were used. The 
crack velocity and the material density were kept constant. 

Figure 8 and Figure 9 show work and internal energy 
variation, respectively. The results reflect higher work and 
higher internal energy for a lower elastic modulus ratio. A 
softer material which has a lower elastic modulus ratio is 
subjected to a larger displacement under the same load while 
all other parameters held the same. This produces higher work 
and higher internal energy since they are proportional to 
displacements. As will be seen in all cases, the trends of 
work and internal energy are always similar to each other. 

Figure 10 shows the kinetic energy variation. Kinetic 
energy is higher for a lower elastic modulus ratio, because 
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for this material the material velocity is higher due to the 
larger displacement. Figure 11 shows the fracture energy 
variation. Fracture energy is also higher for a lower elastic 
modulus ratio. 

The increments of work and energies between the initial 
time and a later time are plotted in Figure 12. The larger 
fracture energy for a lower elastic modulus ratio was caused 
by a much larger increase of work relative to the increases of 
internal and kinetic energies for this material. 

In the second case the crack velocity was varied. Three 
different crack velocities were used. The elastic modulus 
ratio and the material density were kept constant. Figure 13 
and Figure 14 show work and internal energy variations, 
respectively. The time needed for a crack propagation to the 
same crack size increment is longer for a lower crack 
velocity. With a longer time, the crack propagation influences 
the rest of the body more widely. The wider influence results 
in a larger displacement. As a result, it produces higher work 
and higher internal energy. A typical displacement variation 
at a point, shown in Figure 15, clearly is larger for a lower 
crack velocity. 

Figure 16 shows the kinetic energy variation. Kinetic 
energy is higher for a higher crack velocity. The reason is, 
as shown in Figure 17, a higher crack velocity produces higher 
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speeds in the material. Figure 18 shows the fracture energy 
variation. Fracture energy is higher for a lower crack 
velocity because a lower crack velocity caused a larger 
increase of work and a less increase of kinetic energy than a 
higher crack velocity. 

The material density was varied as the last study. Three 
different material densities were used. The elastic modulus 
ratio and the crack velocity were kept constant. Figure 13 and 
Figure 14 show the work and internal energy variations, 
respectively. For the same magnitude of an external force, a 
smaller density causes a higher velocity and larger 
displacement of the plate. Because the displacement is larger 
for a lower material density, work and internal energy are 
higher. Figure 22 confirms it. Figure 23 also shows a higher 
speed for a lower material density. 

Figure 24 shows the kinetic energy variation. Kinetic 
energy is higher for a higher material density. The velocity 
is higher for a lower density, but the mass is larger for a 
higher density. Because the kinetic energy is proportional to 
the mass as well as the square of velocity, which case gives 
a higher kinetic energy depends on which term is dominating. 
In this situation large mass yields a higher kinetic energy 
although the material velocity is lower. Figure 25 shows the 
fracture energy variation. Fracture energy is higher for a 
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lower material density because of the same reason as that for 
the case of different crack velocities. The increments of work 
and energies are also plotted in Figure 26. 
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Figure 7 Work and Energy Variation During a Crack Propagation 
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Figure 8 Work Variation for Different Elastic Modulus Ratios 
and Fixed Crack Velocity and Material Density 
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Figure 9 Internal Energy Variation for Different Elastic 
Modulus Ratios and Fixed Crack Velocity and Material Density 
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Figure 10 Kinetic Energy Variation for Different Elastic 
Modulus Ratios and Fixed Crack Velocity and Material Density 
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Figure 11 Fracture Energy Variation for Different Elastic 
Modulus Ratios and Fixed Crack Velocity and Material Density 
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Figure 12 Work and Energy Differences Between the Initial Time 
and a Late Time for Different Elastic Modulus Ratios and Fixed 
Crack Velocity and Material Density 
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Figure 13 Work Variation for Different Crack Velocities and 
Fixed Elastic Modulus Ratio and Material Density 
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Figure 14 Internal Energy Variation for Different Crack 
Velocities and Fixed Elastic Modulus Ratio and Material 
Density 
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Figure 15 Displacement Values for Different Crack Velocities 
and Fixed Elastic Modulus Ratio and Material Density 
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Figure 16 Kinetic Energy Variation for Different Crack 
Velocities and Fixed Elastic Modulus Ratio and Material 
Density 
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Figure 17 Speed Values for Different Crack Velocities and 
Fixed Elastic Modulus Ratio and Material Density 
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Figure 18 Fracture Energy Variation for Different Crack 
Velocities and Fixed Elastic Modulus Ratio and Material 
Density 
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Figure 19 Work and Energy Differences Between the Initial Time 
and a Late Time for Different Crack Velocities and Fixed 
Elastic Modulus Ratio and Material Density 
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Figure 20 Work Variation for Different Material Densities and 
Fixed Elastic Modulus Ratio and Crack Velocity 
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Figure 21 Internal Energy Variation for Different Material 
Densities and Fixed Elastic Modulus Ratio and Crack Velocity 
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Figure 22 Displacement Values for Different Material Densities 
and Fixed Elastic Modulus Ratio and Crack Velocity 
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Figure 23 Speed Values for Different Material Densities and 
Fixed Elastic Modulus Ratio and Crack Velocity 
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Figure 24 Kinetic Energy Variation for Different Material 
Densities and Fixed Elastic Modulus Ratio and Crack Velocity 
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Figure 25 Fracture Energy Variation for Different Material 
Densities and Fixed Elastic Modulus Ratio and Crack Velocity 
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Figure 26 Work and Energy Differences Between the Initial Time 
and a Late Time for Different Material Densities and Fixed 
Elastic Modulus Ratio and Crack Velocity 
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IV. 



CONCLUSIONS AND RECOMMENDATIONS 



The developed computer code has been proven to work well 
for dynamic crack problems. Very good agreements were obtained 
for the known solutions. The parametric study of dynamic crack 
propagation in composite plates provided the following 
conclusions . 

A softer composite plate which has a lower elastic modulus 
ratio showed a higher fracture energy if the crack velocity, 
the material density and other characteristics of the problem 
were held constant. A composite plate with a lower crack 
velocity and a composite plate with a lower density showed a 
higher fracture energy provided all other parameters were held 
constant, respectively. 

In this study a discontinuous crack release technique was 
used to model the crack propagation. This causes sudden shock 
waves which in turn result in oscillations in solutions. A 
technique that releases the propagating crack nodes 
continuously, for example as shown in Ref. 12, is recommended 
to achieve more stable solutions. 
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APPENDIX 



The general stiffness matrix using a bilinear rectangular 
element and the general case of material property [D] matrix, 
that is. 



[D] 



Al 


a 2 


Aa 


An 


A 2 


As 


Ai 


a 2 


Aa. 



is given as below: 



*ii 


*12 


*13 


*14 


*15 


*16 


*17 


*18 


*21 


*22 


*23 


*24 


*25 


*26 


*27 


*28 


*31 


*32 


*33 


*34 


*35 


*36 


*37 


*38 


*41 


*42 


*43 


*44 


*45 


*46 


*47 


*48 


*51 


*52 


*53 


*54 


*55 


*56 


*57 


*58 


*61 


*62 


*63 


*64 


As 


A 6 


*67 


As 


Ai 


*72 


*73 


*74 


As 


Ae 


a 7 


As 


*81 


*82 


*83 


*84 


*85 


*86 


*87 


*88. 



where 



Kii = — D lx bc 3 +& (D 31 +D 13 ) b 2 c 2 + — D 33 b 2 c 
K 12 = D 13 bc 2 +4 (D 12 +D 33 ) b 2 c 2 + -^-D 32 b 2 c 
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^3 = -^^ii^ 3 +4(- J D 31+ D 1 3) J b 2 C 2 + |D33 J b 3 C 



K 14 = -^-D 13 bc 2 +4 (D 12 -D 33 ) b 2 C 2 + jn 32 b 3 C 



K 15 =-jD 11 bc 3 +4 (-D 31 -D 13 ) b 2 C 2 -jD 33 jb 3 C 



K 16 = -^D 13 bc 3 +4(-D 12 -D 33 )b 2 C 2 ~^D 32 b 3 C 



^i7=|^ii^ 3+ 4 (D 31 -D 13 ) b 2 C 2 -^-D 33 b 3 C 



K ia = ^D 13 bc 3 +4(-D 12 +D 33 )b 2 C 2 -lfD 32 b 3 C 



*21 = -y- D 31 bc 3 +4 ( D 21 +D 33 ) b- 2 C 2 + M D 23J b 3 C 



*22 = ^^ 33 ^ 3 + 4 (^ 32 + ^ 23)^ 2 ^ 2 + ^ J -^ 3 ^ 
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K 22 =-^-D 21 bc 3 + 4 ( ~d 21 +d 22 ) b 2 c 2 + ^D 22 b 3 c 
K 2 t = -lyD 22 bc*+UD 22 -D 22 )b 2 C 2 + ^D 22 b*C 

K 25 = ~^D 21 bc 3 +4 (~D 21 ~D 22 ) b 2 c 2 -^D 22 b 3 C 



^26 = -|^ 33^ 3+4 (-D 22 -D 22 ) b 2 C 2 ~^D 22 b 3 C 



^27 = |^3i ^ 3+4 (^21-^33) b 2 c 2 -^-D 22 b 2 c 



K 2a = ^D 22 bc 2 ^ (~D 22 + D 22 ) b 2 c 2 ~^D 22 b 2 c 



K 2X = -^-D 21 bc 2 + 4 ( D 21 -D 12 )b 2 C 2 + jD 22 b 3 C 



K 32 = ^j-D 13 bc 3 +4 ( ~D 12 +D 22 ) b 2 C 2 + -^D 22 b 3 C 
K 2 3 = -y- D 12 bc 3 +4 ( -D 2 2 -D 12 ) b 2 c 2 + -y- D 3 3 b 3 C 
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A7 34 — D 12 bc 2 + 4( D 12 D 22 ) b 2 c 2 + — D 22 b 3 c 



K 2 5 = -| D^iDC 3 +4 ( -Dj 1 +D 13 ) b 2 c 2 - M d 33 Jb 3 c 



^36 = |^13^ 3+4 (^ 12 -^ 33 ) b 2 C 2 -^-D 22 b 2 C 



K 21 = -^D 21 bc 2 +± (D 21 +D 12 ) b 2 c 2 -^D 22 b 2 c 



^38 = -|^i3^ 3+4 (D 12 +D 22 )b 2 C 2 -^D 22 b 2 C 



K 41 = ~lyD 21 bc 2 +4 (D 21 -D 22 ) b 2 C 2 + ^D 22 b 2 C 



K i2 = -lj-D 22 bc 2 +±{-D 22 +D 22 )b 2 C 2 + ^D 22 b 2 C 



K i2 = ^-D 21 bc 2 ^ {-D 21 -D 22 ) b 2 C 2 ^D 22 b 2 C 
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^44 -^3 3 -fc* C 3 + 4 ( ^32~-^23^ c2 + ^22^^ C 



*4 5 = ' -|- D 3 1*> C 3 +4 ( "^2 1 +D 33 ) b 2 C 2 - D 22 b 3 . C 



K^ = ^D 22 bc 2 ^ (D 32 ~D 22 ) b 2 C 2 -^-D 22 b 2 C 



K i7 =-^D 31 bc 2 +4 (D 21 +D 22 ) b 2 C 2 -^D 22 b 2 C 



K< 9 = -^D 32 bc 2 +4 (D 22+ D 22 ) b 2 C 2 -^D 22 b 2 C 



K sl = -^D 11 bc 2 +i(-D 31 -D 12 )b 2 c 2 -^D 23 b 2 c 



K 52 = -^D 13 bc 2 +4(-D 12 -D 33 )b 2 C 2 -^D 32 b 2 C 



K sz = \ D n bc 2 +4 (D 31 -D 13 )b 2 c 2 -^j-D 33 b 2 c 
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*54 = ' f D lZ bc 3 +4 ( ~ D 12 +D 32 ) b *' <=' 2 - -y-- D 22 b ^ C 



JC 5 5 = ^-Z?ii J bC 3 + 4(D 3 i + i? 1 3)i 5 2 C 2 + M D 33 J b 3 c 



K 56 =^-D 13 bc 2 + 4 (D 12+ D 33 ) b 2 C 2 *^-D 32 b 2 C 



K.^-^-D^bc 2 ^ (~D 31+ D 13 ) b 2 C 2 + ^-D 33 b 2 c 



K s& = -lfD 13 bc 2 + 4 (D 12 -D 33 ) b 2 c 2 + ^D 32 b 2 C 



K 61 = -^D 31 bc 2 +4(-D 21 -D 33 )b 2 C 2 -^D 23 b 2 c 



K 62 =-^D 33 bc 2 +4 (-D 32 -D 23 ) b 2 c 2 -^D 22 b 2 c 



K 63 =^D 31 bc 2 ^(D 21 -D 33 )b 2 C 2 ~^-D 23 b 2 C 
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*64 = | D 33^ 3 + 4 ( -D 22 +D 22 ) b 2 C : 2 - M D 22 £ 3 C 



*6 5 = ^ D 31^ 3+4 ( Z:) 2 l + - D 33) jb2 ^ 2+ -y-^ 2 3 jb3c 
*66 = ^ D 33^ 3+4 (°3 2 + ° 2 3)^ 2 ^ 2 + -y-^ 3 ^ 



*67 = - -y- ^31^ 3 +4 ( -^1 +D 3 3 ) ^ 2 ^ 2 + ^ ^ 3 C 



*6 8 = - y - D 33^ 3 +4 ( ^3 2 "^3 ) ^ 2 ^ 2 + 1 ^ 3 C 



*7i = |o ilJ bc 3 + 4(-D 31+ D 1 3) J b 2 c 2 -^-D33 J b 3 c 



*7 2 = |^13^ 3+4 (^l 2 -^33)^ 2 ^ 2 -^^3 2 i> 3 ^ 



* 73 = -^ llJ bc 3 +4 (D 31 +D 13 ) J b 2 c 2 -|D 3 3J b 3 c 
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K,^-^D 13 bc 2 ^{D 12 ^D 33 )b 2 c 2 -^D 32 b 2 c 
K JS = -^j-D 11 bc 2 +4: (D 31 -D 13 ) b 2 c 2 + ^D 33 b 2 c 
K? 6 = -^f D i3bc 2 +4 (-Oi 2 +D 33 ) b 2 c 2 + ^D 32 b 2 c 
Kn^-j-D^bc 2 +4 (-D 31 -D 13 ) b 2 c 2 + -^D 33 b 2 c 



Kjs = lfDi3bc 2 +4(-D 12 -D 23 )b 2 c 2 + ^D 32 b 2 c 



K sl = -^D 31 bc 2 +4(-D 21 +D 33 )b 2 C 2 -lj-D 23 b 2 C 



K S2 = ^D 33 bc 2 +4(D 32 -D 23 )b 2 c 2 ~lfD 22 b 2 C 
Ks 3 = -^D 31 bc 2 +4(D 21 +D 33 )b 2 C 2 -^D 23 b 2 C 
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K 84 = --^ D 33 bc2+A ( D 32 +D 2 3 ) b 2 C 2 ~^D 22 b 2 C 



*85 = - 



16 



D 3lJ bc 3 +4(D 2 ,-D. 



33 



) b 2 c 



2 +—d. 



23 



b 2 c 



K 86 = -^fD 23 bc 2 +U-D 22 +D 23 )b 2 C 2 + £D 22 b 2 C 



K sl =^-D 21 bc 2 +±(-D 21 -D 22 )b 2 C 2 + ^D 22 b 2 C 



K 8S = ^f-D 2:i bc 2 +4(-D 22 -D 23 )b 2 C 2 + ^D„b 2 C 



' 22 ■* 



In this study the [D] matrix for a unidirectional composite 
material for the plane stress conditions is; 



[D] 





0 

0 
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and the nonzero elements of the [D] matrix are; 






0i 

1- U 12 U 21 



022 ' 



1- U 12 U 21 



„ _ U 1202 _ U 2101 
U \2 ~21 



^ U 12 U 21 1 U 12 U 21 



033 -G 12 



where subscripts 1 and 2 on the right hand side of the 
equations are the direction of the fibers and the transverse 
direction, respectively. Substituting these values into the 
general stiffness matrix gives the stiffness matrix for a 
unidirectional composite. 
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